This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.
library(sf)
library(here)
st_layers("~/Documents/_UCL_grad school/UCL/Term 1 /GIS CASA/casa_0005_coursework/week 3/gadm36_AUS_gpkg/gadm36_AUS.gpkg")
Driver: GPKG
Available layers:
Ausoutline <- st_read("~/Documents/_UCL_grad school/UCL/Term 1 /GIS CASA/casa_0005_coursework/week 3/gadm36_AUS_gpkg/gadm36_AUS.gpkg",
layer='gadm36_AUS_0')
Reading layer `gadm36_AUS_0' from data source
`/Users/christine/Documents/_UCL_grad school/UCL/Term 1 /GIS CASA/casa_0005_coursework/week 3/gadm36_AUS_gpkg/gadm36_AUS.gpkg'
using driver `GPKG'
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
Geodetic CRS: WGS 84
print(Ausoutline)
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
Geodetic CRS: WGS 84
GID_0 NAME_0 geom
1 AUS Australia MULTIPOLYGON (((158.6928 -5...
st_crs(Ausoutline)$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"
library(raster)
library(terra)
jan<-terra::rast("~/Documents/_UCL_grad school/UCL/Term 1 /GIS CASA/casa_0005_coursework/week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_01.tif")
# have a look at the raster layer jan
jan
class : SpatRaster
dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : wc2.1_5m_tavg_01.tif
name : wc2.1_5m_tavg_01
min value : -46.042
max value : 34.065
#average temperature january (01 month)
plot(jan)
newproj<-"ESRI:54009"
# get the jan raster and give it the new proj4
pr1 <- jan %>%
terra::project(., newproj)
plot(pr1)
#ok that was just testing projections
pr1 <- pr1 %>%
terra::project(., "EPSG:4326")
plot(pr1)
# look in our folder, find the files that end with .tif and
library(fs)
dir_info("week 3/")
library(tidyverse)
listfiles<-dir_info("week 3/wc2.1_5m_tavg/") %>%
filter(str_detect(path, ".tif")) %>%
dplyr::select(path)%>%
pull()
#have a look at the file names
listfiles
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_01.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_02.tif
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_03.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_04.tif
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_05.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_06.tif
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_07.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_08.tif
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_09.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_10.tif
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_11.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_12.tif
worldclimtemp <- listfiles %>%
terra::rast()
#have a look at the raster stack
worldclimtemp
class : SpatRaster
dimensions : 2160, 4320, 12 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
sources : wc2.1_5m_tavg_01.tif
wc2.1_5m_tavg_02.tif
wc2.1_5m_tavg_03.tif
... and 9 more source(s)
names : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ...
min values : -46.042, -44.800, -57.986, -64.200, -64.829, -64.395, ...
max values : 34.065, 32.908, 33.081, 34.277, 36.299, 38.458, ...
worldclimtemp$Jan
class : SpatRaster
dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : wc2.1_5m_tavg_01.tif
name : Jan
min value : -46.042
max value : 34.065
#create data frame of cities in australia
site <- c("Brisbane", "Melbourne", "Perth", "Sydney", "Broome", "Darwin", "Orange",
"Bunbury", "Cairns", "Adelaide", "Gold Coast", "Canberra", "Newcastle",
"Wollongong", "Logan City" )
lon <- c(153.03, 144.96, 115.86, 151.21, 122.23, 130.84, 149.10, 115.64, 145.77,
138.6, 153.43, 149.13, 151.78, 150.89, 153.12)
lat <- c(-27.47, -37.91, -31.95, -33.87, 17.96, -12.46, -33.28, -33.33, -16.92,
-34.93, -28, -35.28, -32.93, -34.42, -27.64)
#Put all of this inforamtion into one list
samples <- data.frame(site, lon, lat, row.names="site")
# Extract the data from the Rasterstack for all points
AUcitytemp<- terra::extract(worldclimtemp, samples)
#create a subset for perth
library(tidyverse)
#define where you want the breaks in the historgram
userbreak<-c(8,10,12,14,16,18,20,22,24,26)
# remove the ID and site columns
Perthtemp <- Aucitytemp2 %>%
filter(site=="Perth")
t<-Perthtemp %>%
dplyr::select(Jan:Dec)
hist((as.numeric(t)),
breaks=userbreak,
col="red",
main="Histogram of Perth Temperature",
xlab="Temperature",
ylab="Frequency")
histinfo <- as.numeric(t) %>%
as.numeric()%>%
hist(.)
histinfo
$breaks
[1] 12 14 16 18 20 22 24 26
$counts
[1] 2 2 2 1 1 2 2
$density
[1] 0.08333333 0.08333333 0.08333333 0.04166667 0.04166667 0.08333333 0.08333333
$mids
[1] 13 15 17 19 21 23 25
$xname
[1] "."
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
AusoutSIMPLE <- Ausoutline %>%
st_simplify(., dTolerance = 1000) %>%
st_geometry()%>%
plot()
print(Ausoutline)
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
Geodetic CRS: WGS 84
GID_0 NAME_0 geom
1 AUS Australia MULTIPOLYGON (((158.6928 -5...
crs(worldclimtemp)
[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
Austemp <- Ausoutline %>%
# now crop our temp data to the extent
terra::crop(worldclimtemp,.)
# plot the output
plot(Austemp)
exactAus<-terra::mask(Austemp, Ausoutline)
#subset using the known location of the raster
hist(exactAus[[3]], col="red", main ="March temperature")
exactAusdf <- exactAus %>%
as.data.frame()
library(ggplot2)
# set up the basic histogram
gghist <- ggplot(exactAusdf,
aes(x=Mar)) +
geom_histogram(color="black",
fill="white")+
labs(title="Ggplot2 histogram of Australian March temperatures",
x="Temperature",
y="Frequency")
# add a vertical line to the hisogram showing mean tempearture
gghist + geom_vline(aes(xintercept=mean(Mar,
na.rm=TRUE)),
color="blue",
linetype="dashed",
size=1)+
theme(plot.title = element_text(hjust = 0.5))
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.
squishdata<-exactAusdf%>%
pivot_longer(
cols = 1:12,
names_to = "Month",
values_to = "Temp"
)
twomonths <- squishdata %>%
# | = OR
filter(., Month=="Jan" | Month=="Jun")
meantwomonths <- twomonths %>%
group_by(Month) %>%
summarise(mean=mean(Temp, na.rm=TRUE))
meantwomonths
ggplot(twomonths, aes(x=Temp, color=Month, fill=Month)) +
geom_histogram(position="identity", alpha=0.5)+
geom_vline(data=meantwomonths,
aes(xintercept=mean,
color=Month),
linetype="dashed")+
labs(title="Ggplot2 histogram of Australian Jan and Jun
temperatures",
x="Temperature",
y="Frequency")+
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))
data_complete_cases <- squishdata %>%
drop_na()%>%
mutate(Month = factor(Month, levels = c("Jan","Feb","Mar",
"Apr","May","Jun",
"Jul","Aug","Sep",
"Oct","Nov","Dec")))
# Plot faceted histogram
ggplot(data_complete_cases, aes(x=Temp, na.rm=TRUE))+
geom_histogram(color="black", binwidth = 5)+
labs(title="Ggplot2 faceted histogram of Australian temperatures",
x="Temperature",
y="Frequency")+
facet_grid(Month ~ .)+
theme(plot.title = element_text(hjust = 0.5))
library(plotly)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attaching package: ‘plotly’
The following object is masked from ‘package:raster’:
select
The following object is masked from ‘package:ggplot2’:
last_plot
The following object is masked from ‘package:stats’:
filter
The following object is masked from ‘package:graphics’:
layout
# split the data for plotly based on month
jan <- squishdata %>%
drop_na() %>%
filter(., Month=="Jan")
jun <- squishdata %>%
drop_na() %>%
filter(., Month=="Jun")
# give axis titles
x <- list (title = "Temperature")
y <- list (title = "Frequency")
# set the bin width
xbinsno<-list(start=0, end=40, size = 2.5)
# plot the histogram calling all the variables we just set
ihist<-plot_ly(alpha = 0.6) %>%
add_histogram(x = jan$Temp,
xbins=xbinsno, name="January") %>%
add_histogram(x = jun$Temp,
xbins=xbinsno, name="June") %>%
layout(barmode = "overlay", xaxis=x, yaxis=y)
ihist
# mean per month
meanofall <- squishdata %>%
group_by(Month) %>%
summarise(mean = mean(Temp, na.rm=TRUE))
# print the top 1
head(meanofall, n=1)
# standard deviation per month
sdofall <- squishdata %>%
group_by(Month) %>%
summarize(sd = sd(Temp, na.rm=TRUE))
# maximum per month
maxofall <- squishdata %>%
group_by(Month) %>%
summarize(max = max(Temp, na.rm=TRUE))
# minimum per month
minofall <- squishdata %>%
group_by(Month) %>%
summarize(min = min(Temp, na.rm=TRUE))
# Interquartlie range per month
IQRofall <- squishdata %>%
group_by(Month) %>%
summarize(IQR = IQR(Temp, na.rm=TRUE))
# perhaps you want to store multiple outputs in one list..
lotsofstats <- squishdata %>%
group_by(Month) %>%
summarize(IQR = IQR(Temp, na.rm=TRUE),
max=max(Temp, na.rm=T))
# or you want to know the mean (or some other stat)
#for the whole year as opposed to each month...
meanwholeyear=squishdata %>%
summarize(meanyear = mean(Temp, na.rm=TRUE))